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' ^h' a random matrix approach is used to analyze the vibrational properties of amor- 

C phous solids. We investigated a dynamical matrix M = AA^ with non-negative 

eigenvalues e = oJ^. The matrix A is an arbitrary real N x N sparse random matrix 
with n independent non-zero elements in each row. The average values {Aij) = 
and dispersion (^^j^ — for all non-zero elements. The density of vibrational 
^ states gioj) of the matrix M for A, n ^ 1 is given by the Wigner quarter circle law, 

g{oj) = {N / TmV'^)\/ AnV^ — w^, with radius independent of N. We argue that for 
Q <^ N this model can be used to describe the interaction of atoms in amorphous 

ifi solids. The level statistics of matrix M is well described by the Wigner surmise and 

^ corresponds to repulsion of eigenfrequencies. The participation ratio for the major 

^ part of vibrational modes in three dimensional system is about 0.2 — 0.3 and inde- 

00 pendent of N. Together with term repulsion it indicates clearly to the delocalization 

I of vibrational excitations. We show that these vibrations spread in space by means 

of diffusion. In this respect they are similar to diffusons introduced by Allen, Feld- 
Q man, et al., Phil. Mag. B 79, 1715 (1999) in amorphous silicon. Our results are 

^-H in a qualitative and sometimes in a quantitative agreement with molecular dynamic 

J> simulations of real and model glasses. 

X 
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I. INTRODUCTION 

One of the main problem in the physics of disordered systems is estabhshing the 
general properties of vibrations in amorphous solids (glasses). At low frequencies u there 
are usual plane wave excitations (longitudinal and transverse phonons) characterized by 
well defined wave vector k. These are delocalized excitations. Their density of states 
(in three dimensional space) gph{uj) oc u"^. 

According to recent investigations besides phonons in glasses exist low frequency 
quasilocal vibrational modes^. Their density of states, gqw{u}), is a universal function 
of u. At low frequencies, below some characteristic frequency Ub, the density of states 
5'qiv(i^) oc cj^, and at higher frequencies for u > Ub, 5'qiv(w) oc u. As a result the reduced 
density of states gc^\v{uJ)/u'^ has a broad maximum at frequency u ^ Ub- This is the 
boson peak which was found almost in all glasses. Its theory was developed recently 

iiPa. 

However the range of the boson peak in glasses normally occupies only 10% of vibra- 
tional spectrum. At higher frequencies, for example in amorphous Si02, the vibrational 
density of states g{u) approaches a constant value which extends towards maximum 
frequencies (of the order of the Debye frequency lodJ^. The same behavior of g{uj) was 
found in the soft-sphere glas^, and in amorphous S^. In other glasses the density of 
states has a broad maximum around 00 d and then decays to zercP^. 

The nature of these higher frequency excitations which dominate in the vibrational 
spectra is not yet known and has been discussed a lot in the literatur d^° * ^^ l The main 
question is whether these vibrations are plane waves (phonons) or not? And if they 
are not phonons, whether they are localized or delocalized? The solution of this cru- 
cial question is very important, for example, for elucidating the nature of the thermal 
conductivity of amorphous dielectrics in this frequency (or temperature) range. 

For example the numerical studies of amorphous silicon!^ show that the lowest 4% 
of vibrational modes are plane-wave like {propagons, or phonons) and the highest 3% of 
modes are localized (locons). The rest are neither plane- wave like nor localized. They 
were called diffusons, having in mind that they spread in space by means of diffusion 
(from one atom to another)!^. Many computer simulations of real and model glasses 
show that the most part of vibrational modes are indeed delocalizecP'^. 

One of the goals of this work is to answer the crucial question what is the nature 
of vibrational excitations in disordered systems and in particular of diffusons using the 
methods of random matrix theory. These methods were very successful in understanding 
of universal properties of electronic disordered systems, in theory of financial markets, 
random networks, wireless communications, etc. However there are very small amount 
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of papers using these methods for investigation of vibrations in disordered systems (see 
for exampl^I^. 

We would hke to fill this gap. However we do not want to attach ourselves to some 
specific disordered system (amorphous solid, liquid or glass). We will consider the dy- 
namical vibrational matrix of a general form which is a symmetric random matrix with 
positive eigenvalues. Therefore we hope that properties of this matrix will be sufficiently 
general and independent of particular details of vibrational system. Some of these prop- 
erties may resemble those in real glasses. 

In the present paper we restrict ourselves to vibrational disordered systems where 
are no delocalized low frequency plane wave excitations (phonons). Therefore if vibra- 
tional modes in our system will happen to be delocalized it would have no relation with 
phonons. In fact for existence of phonons it should exist some order in the disordered 
system. Therefore, we consider this case separately somewher^^. 



II. RANDOM MATRIX APPROACH TO VIBRATIONS IN DISORDERED 

SYSTEMS 

Newton equations of mechanical vibrational system in harmonic approximation can 
be presented in a form (scalar modeV^ 



N 

z = l,2,...,iV. (2.1) 

i=i 



Here rrii and Ui are masses and atomic displacements correspondingly, $jj is a force con- 
stant matrix (Hessian). It is a real, symmetric and positive definite matrix. Introducing 
the new variables 

Qi = Uiy/m~i, Mij = ^ij/y/niirrij (2.2) 

we come to the equations 

TV 

w'gi = ^M,,g,-, 2 = 1,2,..., AT. (2.3) 

i=i 

The dynamical matrix M, as well as matrix $ is a real, symmetric and positive definite 
matrix. At some values e = u"^ which are eigenvalues of the djTiamical matrix, the 



system (2.3) has non-zero solutions. These values of u are vibrational eigenfrequencies 
of the mechanical system. 

To describe vibrations in amorphous solids we are going to apply methods of random 
matrix theory. At first glance it is a rather difficult mathematical problem. Not every 
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random symmetric dynamical matrix is appropriate for studying vibrations. Matrix M 
should be also positive definite. It ensures mechanical stability of the system. For that, 
matrix elements Mij should be correlated with each other in a non trivial way. Therefore 
they cannot be taken as independent random numbers. We will solve this problem using 
a following mathematical approach. 

Every real symmetric and positive definite matrix M one can always present in the 
formP^. 

M = AA^. (2.4) 

Here A is some real (not necessary symmetric) matrix of a general form. And vice versa 
for every real matrix A the product AA^ is a positive definite symmetric matrix!^. One 
may assume that in amorphous solids (in some frequency interval) e = uj'^ are eigenvalues 
of matrix AA^ where A is taken in some way randomly. 

Distribution of eigenvalues for matrix M = AA^ for some classes of random matrices 
A was first derived in the paper of Marchenko and Pastur^^. Then it was investigated 
in the theory of financial market^, complex network^ and wireless communication^. 
As far as we know for vibrations in disordered solids this approach was not used so far 
(as an exception see papei^. In the paper^^^'^ it was mainly investigated the case 
of Wishart ensemble^^' when matrix A is a random matrix with independent matrix 
elements with zero mean {Aij) = and equal dispersion (^4^^) = V^. 

It was shown in these papers that for a real matrix A (dimensions x A^) the 
distribution of frequencies u (square roots from eigenvalues e of the matrix M = AA'^) 
for A^ ^ 1 has a quarter circle form 

g(co) = -^V4:NV^ - 0<co< 2VVn. (2.5) 

It formally coincides with the well known Wigner semicircle^' for distribution g{e) of 
the eigenvalues e of symmetric random matrix H with independent random matrix 
elements with zero mean (Hij) = and equal dispersion {Hfj'^ = V^. In this case 
Eq. (2.5) with additional coefficient 1/2 and replacement w to e is valid in the interval 
-2VVN <e < 2Vy/N. 

In the random medium model described by the Wishart ensemble, each of the elements 
Mij of the dynamical matrix M is not zero 

M,j = ^AikAjk. (2.6) 

k 

Obviously this corresponds to the case of long-range interaction, when each atom is 
connected by random springs with every atom in the system. However, this model is not 
justified from the physical point of view. In amorphous materials, only closely spaced 
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atoms are bound together by elastic springs. Therefore the more real would be the 
case when the number of nonzero elements m in each row (and column) of the matrix 
M is small as compared to and does not depend on A^. As a result, the matrix 
M should be sparse. Actually, such sparse matrices arise in computer calculations of 
atomic vibrations in amorphous solids (and liquids). For example, in the case of the 
short-range order for a simple cubic lattice with the interaction only between nearest 
neighbors and the vector character of vibrations (in three dimensional space), we have 
m = 18 + 3 = 21. For other lattices we have m = 24 + 18 + 3 = 45 for a bcc lattice and 
m = 36 + 18 + 3 = 57 for a fee lattice. In the last two cases we took into account all 
interactions in the first and second coordination shells. 

Therefore we get a more real case if we consider a sparse matrix A where each row 
contains only n nonzero matrix elements (with n ^ A^). Then, each row of the matrix 
M = AA^ will have, on average, m = nonzero elements. At rt^ <^ A^, this corresponds 
to the case of a sparse matrix M. We will show below, that in this case, for n ^ 1, the 
density of states is also described by the quarter-circle law 

. . A^ 



V4:nV^ - u^, 0<u < (2.7) 



but with radius independent of A^. The derivation of formulas (2.5) and (2.7) is given 
below, because in this form (and as applied to the problem of vibrations) it is absent 
in the literature. The derivation given irf^ is mathematically rather difficult for non- 
specialists. 



III. THE QUARTER-CIRCLE LAW 



First let us derive Eq. (2.5). We assume that A is a random square matrix N x N 
(with A^ ^ 1) of a general form (not necessarily symmetric) with independent random 
matrix elements Aij so that for any i and j. 

(A,) = 0, (A,/) = V''. (3.1) 

To derive the distribution density of eigenvalues for matrix M = AA^ we will use the 
method similar to the derivation of the distribution density of eigenvalues of a random 
symmetric matrixP^. 

For that we fix matrix A and add to it an independent small random matrix SA with 
all its elements being independent so that for any i and j 



{6Aj) = 0, {SAl)=v\ 



(3.2) 
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with V <^ y. Then the matrix M is changed by 

SM^{A + SA) ■ {A + SAf - M = 6A- A^ + A-SA^ + SA- SA^. (3.3) 

Let U be an orthogonal matrix that transforms the matrix M into a diagonal form. In 
this representation the new matrix 5M can be defined as 

5M = U-^ -SM -U. (3.4) 

According to the perturbation theory the change of the iih. eigenvalue £j of the matrix 
M (due to the addition of the matrix 5M) can be written in the form 

~ 2 

Let us average this equation over random matrix bA. The matrix A we keep fixed 
therefore C/, M, and eigenvalues Si are also fixed. For this purpose we have to calculate 
the mean value and the variance of the elements of matrix 8M. The mean value of the 
diagonal elements is given by 

IbMi^ = ( ^kiUli[Ak^bAim + Mfc^Am + 5AfemM;m) )• (3.6) 

The first and second terms in parentheses are equal on average to zero since = 0. 
The third term at /c 7^ / is also equal on average to zero but at A; = / it is equal to . 
As a result we get 

Im,^ = J2 Ul {5AuJ) = Uh = v'N. (3.7) 

km km 

In the last equality we used orthogonality of the matrix U. 
The variance of the off-diagonal elements {i 7^ j) is given by 

{^fLW,/^ = ^ ^ UkiUij{AkmSAim + SAkmAim + Mfc^M,,,,)) (3.8) 
klm 

We can neglect the last term in the parentheses because it leads to the higher order 
corrections. Then 

(^SMij^^ ^{(Yl UkiUijAkmSAi^y^ + ( ( UkiUuSAkmAi^y) + 

klm klm 

+2(^Y UkiUljAkmSAim ■ Y UkiUij6AkmAim). (3.9) 
klm klm 
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In the first term we expand the square. The nonzero contribution is given only by the 
terms that contain the product of the same elements of the matrix 6A. Their mean is 
equal to f ^ 

2\ 



^ ^ UkiUljAkfn6Aim 



Mm 



^ ^ Uk-iiUk2iUlj Ak^rnAk2m'V 



(3.10) 



The sum over / of Uij^ gives 1 due to orthogonality conditions for matrix U. The 
reminder is the diagonal element (u) of the matrix U^^AA^U which is the eigenvalue 
Ei- Therefore the first term in Eq. (3.9) is equal to v'^Ei. By renaming the indices we 
similarly find that the second term in Eq. (3.9) is v'^Ej. Because of orthogonality of the 



matrix f/, the third term in Eq. (3.9) is equal to zero. As a result we get 

2\ 



(3.11) 



Let us now insert Eqs. (3.7, 3.11) into the right-hand side of averaged Eq. (3.5). As a 
result the mean value of the correction to the ith eigenvalue can be written in the form 
(we will keep the same notation for it) 

2 X ^ ~\~ 



Se {e,,v')=v'N + v'J2-^ 



Taking into account that ^ 1, we include the first term into the sum 

1 



6e (e,, V^) = 2vhi J2 



(3.12) 



(3.13) 



For A^ ^ 1, the distribution density a{E, V"^) of eigenvalues of the matrix M is a contin- 
uous function of E. Then sum ( 3.13[ ) can be approximated by the principal value of the 
integral 

(3.14) 



bE {e,V^)='2v\ I ^j^dE'. 



Now let us consider a change of the number of eigenvalues of the matrix M in the 
interval {e, e + Ae) due to the change of the matrix A to the matrix A + 6 A. On the one 
hand this change is given by 

d{a6E) 



a {e, 1/2) 6e {e, V^) -a{E + Ae, V'^) 6e {e + Ae, V^) 



dE 



-Ae. 



(3.15) 



On the other hand since we added an independent random matrix 5A to the matrix A, 
the variance of the elements of the new matrix A + 5 A increased from V'^ to V'^ + v'^. 
Therefore the change of the number of eigenvalues in the interval [e, e + Ae) can be 
presented in the form 



a [e, + i;2) AE-a [e, V^) Ae 



2 da ^ 



(3.16) 
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By equating (3.15) and (3.16) we obtain the second equation 

d{a5e) 



de 



> da 



(3.17) 



To solve the system of two equations (3.14) and (3.17) let us introduce the new vari- 
able X = e/V"^. Since eigenvalues of the matrix M are proportional to V"^, the quantity 
X = e/V'^ for each eigenvalue remains unchanged when all elements of the matrix A 
are multiplied by a constant. Secondly, according to the normalization condition, the 
density of eigenvalues is scaled as l/V"^. As a result the distribution density can be 
presented in the form 

^(^'^') = ^^(^)' (3-18) 



where a{x) is some function of x. Making use of it we find from Eq. (3.14) that 6e can 
be presented in the form 

6e{e,V^) =v^6e{x). (3.19) 



Then Eq. (3.14) transforms as follows 



6eix) = 2x I ^^dx'. 



X — X 



(3.20) 



Substitution of Eqs. (3.18, 3.19) into Eq. (3.17) gives us the differential equation with 
one variable 

d \ a5e 

(3.21) 



d[xa) 



dx 



dx 



It follows from Eq. (3.20) that fe(0) = 0. As a result after integration of Eq. (3.21) we 

get 

5e{x) = X. (3.22) 



Finally using that we derive from Eq. ( 3.20[ ) the integral equation 

X — X 2 

The normalization condition can be written in the form 



(3.23) 



(3.24) 



The solution of Eq. (3.23) with the normalization condition (3.24) is given by 
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Coming back to the original variables we get the distribution density of eigenvalues for 
the matrix M 

"(^•''''"SSTj/^^T^- 0<^<''W=- (3-26) 
This distribution is a particular case of the Marchenko-Pastur distributior l^^ * ^^ l Taking 
into account that e = u'^ the density of vibrational states can be represented as 

g(uj) = -^V^NV^ - < w < 2VVn. (3.27) 

This equation has the form of a quarter-circle. 

As was already mentioned above in Section [IT] this model with the long-range inter- 
action is not justified from the physical point of view for vibrations in amorphous solids. 
On the one hand it implies that each atom interacts with every atom in the system. On 
the other hand it gives the width of the vibrational spectrum depending on the size of 
the system N. In the next section we consider the sparse random matrix A. In such a 
way one can overcome these two limitations. 



IV. A SPARSE RANDOM MATRIX 



Let us consider a sparse random matrix A. Now in each row i we randomly choose 
n positions j and write random numbers Ojj there with the same probability density 
pQ{aij). In all other positions we write zeros. As before all nonzero random matrix 
elements are characterized by their mean (aij) = and the variance (afj) = V^. An 
example of such sparse random matrix A (8 x 8) with n = 3 is shown below 

/oo*o*oo*\ 


A= 00*0*0*0 . (4.1) 



\0*00**00/ 

The asterisks here indicate nonzero random matrix elements. 

The elements of such random matrix are independent from each other and their 
distribution is given by the probability density 

p(x) = ^poix) + ^^^Six). (4.2) 

Here S{x) is the Dirac delta function. The mean value of matrix elements (Aij) = 0. 
The variance of this distribution is given by 

^^P{x) d^ = J^ j ^Vo(a:) dx = -V\ (4.3) 
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The averaging is performed over both zero and nonzero elements of the matrix A. 



Such sparse random matrix A satisfies to conditions similar to Eqs. (3.1) used for 



derivation of Eq. (3.27). Therefore, for n ^ 1 we can use the derived quarter-circle 
density of states (3.27) by substituting the variance nV'^/N instead of the variance V^. 
As a result 

N I 

(4.4) 



n. 



At the same time we can take n N. This means that for N ^ n ^ 1 the quarter-circle 
distribution for the density of states g{uj) is still valid even in the case when nonzero 
elements occupy only a small part of the matrix A. Such form of the distribution g{u) 
in our model is a universal law and does not depend on the distribution density Po{aij), 
the size of the system N, and the number of nonzero elements n for sufficiently large 
values of n. 




FIG. 1: Numerically calculated density of states g{u}) of the matrix M = AA^ for two different 
values of n and N = 1000. Black line shows the density of states given by Eq. ( |4.4[ ). The 
random matrix A has a normal Gaussian distribution of nonzero matrix elements po{aij) with 



{o-ij) = and variance (afj) = 1. 



The numerical analysis confirms that with an increase of n the density of states 
g{u) indeed approaches the quarter-circle distribution (for n ^ 1) and in this case 
the inequality n ^ is possible (see Fig. [T]). Already for n of the order of 10 and 
A^ = 1000, we get the density of states that is only slightly different from the quarter- 
circle distribution. Normalized to unity it does not depend on the size of the system 
N. 

For such sparse random matrix A with n nonzero elements in each row, each row 
(and each column) of the symmetric matrix M has on average nonzero elements (for 
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^ A^). Physically this implies that one atom interacts with approximately close 
atoms. The possibility to take a fixed value of n ^ 1 and allowing N ^ oo gives us 
the opportunity to describe the interaction between close atoms in amorphous materials 
what is physically reasonable. 

From Eq. (4.4) we can obtain the distribution p{e) of the eigenvalues e = of the 
matrix M = AA^ 

. , N UnV^ - e , 

This expression has a singularity at e — )■ oc l/A/e). Probably this singularity 

in the distribution of eigenvalues of the dynamical matrix was observed iri^ for the 
amorphous and liquid phases of Si02. But it was attributed (improperly in our opinion) 
to special correlations between diagonal and non-diagonal elements of the dynamical 
matrix. A similar behavior at e — )■ was observed for the distribution of eigenvalues 
irP^ (see Fig. 2 of this paper). 

In conclusion of this section, it should be noted that the singularity in the density of 
states g{oj) at a; — manifests itself for small values of n (see Fig. [T] for n = 5). Similar 
singularity exists also in the density of eigenvalues g{e) (at e — )■ 0) of the symmetric 
sparse random matrix Hij for small values of n. The singularity is suppressed with 
increase of n when the density of eigenvalues g[e) approaches to the form of Wigner 
semicircl^2SH2l_ xhis singularity was first discovered in the density of vibrational states 
of a disordered one-dimensional chain by DysorP^. Therefore it sometimes is called the 
Dyson singularity. It has been believed that this singularity is an indication of strong 
fiuctuations in a random medium and related localization of mode^. We hope to 
consider this problem in more details in a separate work. 



V. THE CUBIC LATTICE WITH RANDOM BONDS 

The symmetric sparse random matrix M = AA^ considered in the previous section 
is topologically equivalent to a tree (closed to itself on the system size). The number 
m = 71^ specifies the order of branching or the coordination number of this tree. However 
a random bonds structure in amorphous solids (glasses) more likely corresponds to the 
short-range order in the atomic arrangement topologically similar to the bond structure 
existing in the corresponding crystals. It is clear that topologically a crystal structure 
differs fundamentally from a tree structure. In a tree structure there are no closed loops 
of bonds which are present in a lattice. Therefore our purpose in this section is to 
generalize obtained results to spatial structures with the topology of a crystal lattice 
rather than a tree. Otherwise our system remains random without any periodicity 
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(except for the topology of the elastic bonds). 

Let us consider how making use of a random matrix A we can build a dynamical 
random matrix M = AA^ with the known topological bond structure but having random 
bond strengths. As an example we consider a simple cubic lattice wxwxw with N = 
atoms. The atoms have integer coordinates (x, y, z) and each coordinate can take on 
values from 1 to w. Let us introduce the integer index i = x + w{y — 1) + ^^(2; — 1). 
Each atom in the lattice is characterized by its unique index i running from 1 to N. 

We construct a random matrix A in the following way. The element A^ is nonzero 
and random if the ith and jth atoms are nearest neighbors or it is the same atom (for 
i = j). All other elements Aij we take equal to zero. We note that for i ^ j the elements 
Aij and Aji are independent random numbers (the matrix A is not symmetric). As 
a result for a simple cubic lattice we have seven nonzero elements in each row and in 
each column of matrix A (except for the rows and columns corresponding to boundary 
atoms). In the dynamical matrix M = AA^ the element Mij will be nonzero if zth atom 
will be nearest or next nearest to the jth atom (or it is the same atom for i = j). Fig. |2] 
shows the atoms interacting with the central (black) atom. All other atoms in the lattice 
interact similarly with their neighbors. 
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FIG. 2: Schematic diagram illustrating the interaction of atoms in a simple cubic lattice. 
Shown are atoms interacting with the central (black) atom. Yellow atoms have on average a 
less rigid bonds with the central atom. 

We calculated numerically the density of vibrational states g{<^) for this cubic lattice. 
The mean value of nonzero elements of the matrix A was taken equal to zero, {Aij) = 
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and the variance was taken to be equal unity, (^ij) = 1 (Gaussian distribution). The 
results of the numerical calculations are shown in Fig. |3| For comparison, the results 
for the density of states g{u}) of the sparse random matrix M = AA^ with n = 6 are 
also presented in Fig. [3j In can be seen from the figure that the density of states for 
the cubic lattice with random bonds almost coincides with the density of states for the 
sparse random matrix with the corresponding value of n. It can be shown that with 
increase of the radius of the interaction between atoms (i.e. if we take into account 
the interaction with next neighbor shels) the shape of the density of states g{<jj) will 
approach a quarter-circle form. 




FIG. 3: Comparison of the density of states for the model of a simple cubic lattice with 
random bonds with the model of sparse random matrix A with parameter n = 6. In both 
cases iV = 1000. 



VI. PARTICIPATION RATIO 

One of the most important problems in physics of disordered systems is the problem of 
modes localization. As is well known from the seminal paper of AndersorP^ a sufficiently 
strong disorder in the system leads to localization of elementary excitations. To estimate 
the inverse strength of mode localization one usually introduces the participation ratio. 
As a rule the participation ratio is defined as 

P{^) = , (6.1) 

1=1 
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where ei{uj) is the ith coordinate of the eigenvector corresponding to the eigenvalue cu^ 
of the dynamical matrix M. In the case of completely localized mode 

|ei| = 1, 62 = 63 = ... = Cat = 0, -P ~ ^ (6-2) 

the participation ratio P{uj) decreases with increase of the system size. In the case of 
completely delocalized mode 

|ei| = |62| = ... = |67v| = P^l, (6.3) 

the participation ratio does not depend on the system size and is of the order of unity. 
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FIG. 4: Numerical calculations of the participation ratio for the sparse random matrix A 
{M = AA"^) for N = 10000 and n = 10. Nonzero elements of the matrix A have a Gaussian 
distribution with zero mean and unit variance. 

We performed numerical calculations of the participation ratio for vibrational exci- 
tations in our model. The results for the dynamical matrix M = AAJ- in the case of 
sparse random matrix A are presented in Fig. |4} As can be seen from this figure the 
participation ratio -P(w) almost for all frequencies excluding highest and lowest ones are 
in the range 0.2 < -P(i^) ^ 0.3. One can show that it is independent of A^. Therefore 
almost all vibrational modes in this range are delocalized. Qualitatively our plot for 
PiijS) coincides well with the results of numerical calculations of the participation ratio 
for amorphous Si02 using molecular dynamics method^ in the wide frequency range 
< w < 120 meV. 

As it could be expected we found that an increase of parameter n leads to a decrease 
in the number of localized modes and to an increase of the number of delocalized modes. 
The participation ratio for all frequencies approaches the limit equal 1/3 (see below). 
Since these delocalized modes are not plane waves (phonons), they according to the 
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terminology proposed by Allen et al^ are referred to as diffusons. We have shown that 



their spacial structure is random and character of spreading in space obeys to diffusion 
law. Their investigation can shed light on the mechanism of thermal conductivity in 
amorphous solids. This problem will be analyzed in more details in a separate work. 

Similar results were obtained for the cubic lattice with random bonds. In this case the 
participation ratio P{uj) ~ 0.2 and is slightly lower than that in the previous case. But 
it does not depend on the system size N as well, what indicates to delocalization of the 
modes. It is interesting to note that in a two-dimensional (square) and one-dimensional 
lattices with random bonds (constructed in a similar way) the participation ratio is one 
order of magnitude smaller than that in the 3d cubic lattice. By analogy with disordered 
electronic system^ this can indicate to the localization of vibrational modes in these 
low- dimensional structures. 

The numerical values of the participation ratio for vibrational modes -P(w) in various 
glasses according to the data obtained by molecular dynamics methods usually are in 
the range 0.2 < -P(i^) ^ 0.6^"^. This is in a good agreement with the results of random 
matrix theory. Indeed let us assume that the eigenvectors Ciiuj) {i = 1,2, A^) of the 
random matrix M = AA^ (which are unit vectors in A^-dimensional space) 



are isotropically distributed in all possible directions. Then the quantity r = ej{u) will 



N 




(6.4) 



1=1 



be distributed according to the Porter-Thomas laW^ 




(6.5) 



As a result we have 




(6.6) 



and according to Eq. (6.1) the participation ratio is 




(6.7) 



In the literatur^ one can find another definition of the participation ratio for the 
vector model (in contrast to the above scalar model) 
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where index i indicates the index number of the atom {i = 1, 2, A^) and the index a 
stands for the Cartesian projection of the displacement of the ith atom onto the a axis 
(a = x,y,z). In this case under the same assumption that the unit vectors eia{uj) are 
isotropically distributed in 3iV-dimensional space we have in analogy with Eqs. ( 6.6[ ) 



Since 

{el + el + eL)') = (e^) + {e%) + (ef,> + 2 ((eL> {e%) + (e^) (e^) + (el)) , 

(6.10) 



then using Eqs. (6.9) we find that the participation ratio ^3(0;) (6.8) is equal 



Psitu) = 3/5 = 0.6. (6.11) 

In this sense the values of P{uj) = 1/3 for the scalar model and Psiu)) = 0.6 for the vector 
model are equivalent to each other from the standpoint of the random matrix theory. 
The value of Psiu)) ~ 0.6 was obtained in numerical calculations of the participation 



ratio for a soft-sphere glass^^l. Finally making use of Eqs. (6.9) it can be shown that 
participation ratios Pq ~ Psi ~ 0.3 calculated numerically by Jin et al-^ for amorphous 
Si02 are also in a good agreement with theoretical values Pq = Psi = 1/3 that follows 
from the formula (18) of this work. Summarizing we can conclude that the participation 
ratio calculated in different papers for different glasses are in a good agreement with 
predictions of the random matrix theory. 



VII. LEVEL STATISTICS 

The level statistics is another powerful criterion that makes it possible to judge about 
the localization or the delocalization of vibrational modes. If modes are localized their 
frequencies are randomly distributed over the frequency axis without any correlation 
with each other according to the Poisson distribution. For quantitative description let 
us introduce the normalized difference between eigenfrequencies 

(7.1) 



{At,}- 

Here Ao; = Wj+i — Ui is the distance between two neighboring frequencies that corre- 
sponds to the frequency u ^ Ui and (Aw) is the mean distance between these frequencies. 
Then for localized modes the distribution function Z{s) can be presented in the form 



Z{s) = exp(— s) 



(7.2) 
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and described by the Poisson distribution. 

When modes are delocahzed the term repulsion effect takes place and for small values 
of s ^ 1 we have Z{s) oc s. In the random matrix theory for the Gaussian orthogonal 
ensemble the level statistics is described well by the Wigner surmis#^'^ 

Zw{s) = ^ s exp (-^s^) . (7.3) 

As was shown by Plerou et al.'^ this formula also adequately describes the Wishart 
ensemble. 

Now let us examine analytically the problem of vibrational terms repulsion in our case. 
We want to find out how the probability of finding two close neighboring eigenfrequencies 
depends on the distance between them if this distance is much smaller than the average 
distance (i.e. s ^ 1). For this purpose we derive the probability of finding two close 
neighboring eigenvalues Ei and Ej of the matrix M = AA^ that are the squares of the 
corresponding eigenfrequencies. We will restrict our analysis to the case of Wishart 
ensemble considered in Section IIIII 



In Eq. (3.13) the term in the sum that contains the difference between these two close 
eigenvalues in the denominator will be considerably larger than the other terms. Let us 
take into account only this term. Then 

Se{e,,V') = ^. (7.4) 

Ei Ej 

Let u = Ei — Ej > he the difference between two neighboring eigenvalues. Then its 
change due to the perturbation 6 A of the matrix A is given by 

A-v'^E 

6u {u, E, V^) = 6e {Ei, V^) - 6e {ej, V^) = — , (7.5) 

where e is the mean value between eigenvalues Ei and Ej. 

Let zi^UjEjV^) be the density of distribution of the difference z/ between two neigh- 
boring eigenvalues. Similar to the derivation of Eq. (3.17) we get 

- ^ = v^^. (7.6) 
The function z(z/, e, V"^) can be presented in a form 

z{u,e,V') = ^z{x,e), (7.7) 
where x = u/V"^. Substitution of Eqs. ( 7.5[ 7.7) into Eq. (7.6) gives 
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For small values of x the right-hand side of this equation can be replaced by zero. In 

this case its solution for the function 'z{x) will be proportional to x. This in turn means 

that for small v the function z{v^e, V'^) oc v i.e. there is a term repulsion effect. 

In the case of the sparse random matrix A we have calculated numerically the level 

statistics for the dynamical matrix M = AA^ . The result (averaged over all frequencies) 

is shown in Fig. [5j It follows from this figure that for small s ^ 1 we indeed have 

Z{s) oc s what is in agreement with the analytical result. It can also be seen from the 

Fig. |5] that Wigner formula (7.3) adequately describes the level statistics in the case 

of sparse matrices. As was mentioned above such statistics corresponds to the case of 

delocalized modes. Therefore we conclude that majority of the vibrational modes in our 

system are delocalized. This is in a good agreement with the data presented in the Fig.|4] 

for the participation ratio P{uj). In conclusion we note that our results agree well with 

the molecular dynamics calculations of the statistics of vibrational levels for amorphous 
clusterpEl. 




FIG. 5: Numerical calculation of the density of distribution of the difference between neigh- 
boring frequencies i.e. square roots of the eigenvalues of the matrix M = AA^ . The matrix 
A is a sparse random matrix with = 1000 and the number of nonzero elements n in each 
row. Nonzero elements of the matrix A have a Gaussian distribution with zero mean and unit 



variance. The solid line indicates the theoretical Wigner dependence (7.3). 
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VIII. THE DISTRIBUTION OF ELEMENTS OF THE DYNAMICAL 

MATRIX 

In the last section we consider the statistics of nonzero matrix elements of the dy- 
namical sparse random matrix M = AA^ and compare it with available computer 
calculations performed for amorphous systems by molecular dynamics methods. First 
we analyze the distribution of nonzero off-diagonal elements of the matrix M. As follows 
from the Eq. (2.6) an off-diagonal matrix element Mij is the "scalar product" of the ith 
and jth rows of the matrix A. In order for the element Mij at i ^ j to be nonzero it is 
necessary that the column positions of at least one pair of nonzero elements in the ith 
and jth rows of the matrix A coincide. For a sparse random matrix A (i.e. for n ^ A^) 
the probability of coincidence of more than one pair of nonzero elements is negligible. 
Therefore the off-diagonal elements of the matrix M for the most part are equal to zero 
and if they are nonzero they are the products of two independent nonzero elements of 
the matrix A. Then the density of distribution of nonzero off-diagonal elements of the 
matrix M is given by 

/ S{x-^v)piOpiv)d^dv= -p(-)p{t)dt. (8.1) 

Here p{x) is the distribution function of nonzero elements of the random matrix A. In 
the case when the non-zero elements of the matrix A are normally distributed with zero 
mean and variance we have 



Substitution of the probability density (8.2) into the integral (8.1) gives the integral 



representation of the zeroth order Macdonald function 

Pm W = -^.M\x\/V). (8.3) 

We have compared this function with numerical results of the papei'^ where the 
dynamical matrix for a disordered system of atoms interacting through the Lennard- 
Jones potential was calculated using the molecular dynamic methods. It can be seen 



from the Fig. [6] that for the variance V"^ = 0.01 our Eq. (8.3) is in a good agreement 
with the numerical result^. 

The diagonal elements of the matrix M are distributed as a sum of n squares of 
elements from one row of the matrix A. When the elements of the matrix A are normally 
distributed with zero mean and variance V"^ we have well known distribution of 
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FIG. 6: Comparison of the distribution (8.3) (solid line) with the results of the numerical 
calculationJ231 (open circles). The variance is V'^ = 0.01. 



diagonal elements of the matrix M 



exp 



X \ 



(8.4) 



VT (n/2) 

As it can be seen from the Fig. [T] this distribution (for the same choice of the vari- 
ance V'^ = 0.01 and n = 14) also agrees quite well with the results of the numerical 
calculations irp31. 



diag / \ 

Pm - 



0.5 - 




FIG. 7: Comparison of the distribution (8.4) of the diagonal elements of the matrix M = AA^ 



for n = 14 and V"^ = 0.01 (solid line) with the results of the numerical calculation^^!] (open 
circles). 
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IX. CONCLUSION 

We have demonstrated that the dynamical sparse random matrices of the general 
form M = AAl- with nonnegative eigenvalues £ = cu^ can be successfully used for the 
description of sufficiently general properties of the vibrational spectra of amorphous 
solids. Compared to the currently used molecular dynamics methods they have an im- 
portant advantage that the construction of the random dynamical matrix corresponding 
to a stable system requires significantly less efforts than numerical calculations for real 
glasses with specific interatomic interaction potentials. In many cases the results are 
found to be quite similar. 

The main our result is that the density of states g{io) of the vibrational system 
described by the sparse random matrix M = AA^ is given by the quarter-circle law. In 
our model this form of the density of vibrational states is a universal law. It depends 
neither on the density of distribution of non-zero matrix elements of the sparse random 
matrix A (provided that their mean is equal to zero and the variance is finite), nor on 
the size of the system iV, nor on the number of nonzero elements n in each row of matrix 
A for n ^ 1. This result also does not depend strongly on the topology of elastic bonds 
between the atoms. It was shown for disordered systems with topology of a tree (Bethe 
lattice) and for a cubic lattice with random bonds. 

The study of the problem of localization of these vibrational modes in the three- 
dimensional system led us to the conclusion that despite a high degree of disorder the 
majority of the modes arc dclocalizcd harmonic excitations. This is evidenced by the 
values of the participation ratio and the statistics of vibrational levels where the term 
repulsion effect clearly manifests itself. Our results are in a good agreement with the 
results obtained for real glasses by molecular dynamics methods. 

Finally by using the concept of sparse random matrices we calculated the distribution 
of matrix elements of the dynamical matrix M . The results of the calculations are in 
a good agreement with the numerical data obtained by molecular dynamics methods 
for some class of amorphous systems. The aforesaid allows us to draw the important 
conclusion that random dynamical matrices of the type M — AA^ reflect some universal 
characteristics of vibrational spectra of amorphous solids. 

In conclusion it should be noted that the density of states in the form of the quarter- 
circle distribution obtained in our work means that g{uj) tends to a constant value 
at a; 0. However usually in disordered systems (in three-dimensional space) such 
behavior of g{uj) can occur only in the absence of the low frequency acoustical phonons. 
These phonons are weakly decayed plane wave excitations propagating in space with 
velocity of sound. Therefore as a rule the density of states of glasses at low frequencies 
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has a phonon gap, g{u) oc — )■ at w — )■ 0. As we already mentioned these phonon 
modes usually contain not more than 10% of all vibrational modes. At higher frequencies 
we have delocalized vibrations of different type, namely diffusons. They occupy the larger 
part of the spectrum. We believe that we can use random matrix theory to describe 
them. 

Relatively recently in the literature there have been appeared disordered systems 
of a new type. In these systems the width of the phonon gap can be decreased or 
even reduced to zero by decreasing the rigidity of the system. In the latter case it 
was found that the density of vibrational modes g{u) indeed goes to a constant value 
at w — i- 0. Such behavior of g{uj) was observed in computer simulations of disordered 
systems of the type of granular media in the vicinity of the jamming transition point 
when decreasing the density of particles interacting through repulsive forces with a finite 
effective radiud^^'^. A similar behavior was also observed for other computer model^. 
Let us mention also the calculations performed by Trachenko et al.'^ who calculated the 
density of vibrational states in amorphous Si02 in the rigid tetrahedron approximation. 
In this approximation the density of states g{u)) also was found to be finite at w — )■ 0. We 
hope that our random matrix model also will make it possible to describe such systems. 
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